Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.
library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed
library(bayesplot) # needed for MCMC diagnostics
library(DHARMa) # model validation
library(ggdist) # partial plots
library(tidybayes) # partial plots
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans
library(ggpubr) # arranging figures
library(mgcv) # GAMThis data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty
rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |>
mutate(dev.temp = as.factor(dev.temp),
replicate = str_sub(sampleID, -1,-1),
population = factor(population)) |>
# number of observations = 5758
filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
sampleID != "60.LCKM.152.30.1" # 5637 - 76 = 5542
) |>
filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes)
group_by(sampleID) |>
mutate(max_value_index = which.max(rate.output),
row_number = row_number(),
max.rate = max(rate.output),
low.rate = mean(rate.output[order(rate.output)[1:10]]),
net.rate = max.rate - low.rate) |>
filter(row_number <= max_value_index) |>
ungroup() |>
mutate(region = factor(region),
dev.temp =factor(dev.temp),
level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
tank >= 200 & tank <= 299 ~ 2,
tank >= 300 & tank <= 399 ~ 3,
TRUE ~ NA_real_)),
female = factor(female))eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figeda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figFirst we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.
Hypothesis test will include:
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
f.model.null <- bf(rate.output ~ 1,
family = gaussian())
model.null <- brm(f.model.null,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
f.model1 <- bf(rate.output ~ 1 + (1| female) + (1| tank),
family=gaussian())
model1 <- brm(f.model1,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
f.model2 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level),
family=gaussian())
model2 <- brm(f.model2,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 13 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
f.model3 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population),
family=gaussian())
model3 <- brm(f.model3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
f.model4 <- bf(rate.output ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order),
family=gaussian())
model4 <- brm(f.model4,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Output of model 'model.null':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -31469.1 78.1
## p_loo 2.5 0.4
## looic 62938.2 156.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.6 99.6
## p_loo 57.6 3.1
## looic 61325.1 199.1
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.7 99.5
## p_loo 57.6 3.0
## looic 61325.4 198.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.7 99.6
## p_loo 58.1 3.2
## looic 61325.5 199.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model4':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.5 99.4
## p_loo 57.5 3.0
## looic 61325.0 198.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model4 0.0 0.0
## model1 -0.1 0.4
## model2 -0.2 0.4
## model3 -0.2 0.4
## model.null -806.6 48.2
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
rate.priors <- prior(normal(430, 0.25), class="Intercept") +
prior(normal(0, 40), class="b")
prior(student_t(3, 0, 195), class = "sigma")f.model1.p1 <- bf(rate.output ~ 1 +
dev.temp*region + t +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p1 <- brm(f.model1.p1,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.987, p-value = 0.656
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
f.model1.p2 <- bf(rate.output ~ 1 +
dev.temp*region + poly(t, 2) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p2 <- brm(f.model1.p2,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98869, p-value = 0.744
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
f.model1.p3 <- bf(rate.output ~ 1 +
dev.temp*region + poly(t, 3) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p3 <- brm(f.model1.p3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98619, p-value = 0.576
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
## Output of model 'model1.p1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -29955.1 94.0
## p_loo 60.1 3.1
## looic 59910.2 188.0
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30526.2 100.6
## p_loo 58.4 3.1
## looic 61052.4 201.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30519.9 100.2
## p_loo 59.1 3.2
## looic 61039.7 200.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.2]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model1.p1 0.0 0.0
## model1.p3 -564.8 31.2
## model1.p2 -571.1 31.4
f.model1.gam <- bf(rate.output ~ 1 +
dev.temp*region + s(t) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.gam <- brm(f.model1.gam,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.gam, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.gam, ndraws=250, summary=FALSE)
model1.gam_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gam_resids) ; testDispersion(model1.gam_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98442, p-value = 0.6
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.